function  polynom_coeff  = compute_policy_function(Y,X)

%Rescaling inputs and output for a more accurate matrix inversion. 
%The estimated polynomial coefficients are then rescaled back;
           
State_var = [ones(size(X(:,1))), log(X(:,1)), log(X(:,2)), log(X(:,1)).*log(X(:,2))];
temp_max = max(State_var,[],1);
X_var_max = repmat(temp_max,[size(State_var,1),1]);
X_var_rescaled = State_var./X_var_max;


try

    temp_coeff = (X_var_rescaled'*X_var_rescaled)^(-1)*(X_var_rescaled'*Y);
    coeff = temp_coeff.*( 1./temp_max )';

catch %This is if the policy function is constant;

    coeff = [nanmean(Y)  0 0 0]';

end


polynom_coeff = coeff;